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Abstract — We  developed  a  novel  method  for  estimating 
a  single  parameter  (rp)  which  summarizes  the  diastolic 
filling  dynamics  of  the  ventricle  that  governs  the  effects 
of  diastolic  filling  time  on  stroke  volume.  The  method 
involves  the  analysis  of  beat-to-beat  fluctuations  in  heart 
rate,  arterial  blood  pressure,  and  left  ventricular  stroke 
volume,  each  of  which  may  be  measured  noninvasively  in 
humans.  In  order  to  evaluate  the  method,  we  constructed 
a  computational  model  of  the  human  cardiovascular  sys¬ 
tem.  We  applied  the  method  to  data  generated  from  the 
computational  model  and  found  close  agreement  between 
the  estimated  and  actual  rp.  We  also  demonstrated  with 
computational  model  examples  that  the  method  could  be 
utilized  to  monitor  increasing  degrees  of  cardiac  tampon¬ 
ade  and  mitral  stenosis,  two  cardiovascular  disease  pro¬ 
cesses  which  prohibit  the  ventricles  from  filling  adequately. 
This  model-based  study  motivates  the  experimental  vali¬ 
dation  of  the  method. 

Keywords —  system  identification,  cardiovascular  model¬ 
ing,  ventricular  diastolic  compliance 

I.  Introduction 

Cardiac  tamponade,  restrictive  cardiomyopathy,  and 
mitral  stenosis  are  examples  of  disease  processes  which 
prohibit  the  ventricles  from  filling  adequately.  In  this 
paper,  we  present  a  method  for  the  quantitative  charac¬ 
terization  of  ventricular  diastolic  filling  dynamics  which 
requires  beat-to-beat  measurements  of  heart  rate  (HR), 
arterial  blood  pressure  (ABP) ,  and  left  ventricular  stroke 
volume  (SV).  These  measurements  may  be  obtained  non¬ 
invasively  in  humans  with,  for  example,  Doppler  ultra¬ 
sound  and  Finapres  techniques.  The  method  may  there¬ 
fore  potentially  provide  a  powerful  means  for  monitor¬ 
ing  patients  who  are  being  treated  for  these  cardiovas¬ 
cular  disease  processes.  The  method  also  complements 
another  algorithm  that  we  previously  developed  which 
requires  the  same  beat-to-beat  measurements  but  char¬ 
acterizes  the  total  peripheral  resistance  barorcflex  [1]. 

Ventricular  diastolic  filling  dynamics  are  highly  com¬ 
plex  and  are  determined  by  a  multitude  of  parameters. 
This  may  be  illustrated  with  Fig.  1  which  includes  a 
sketch  of  right  ventricular  stressed  volume  as  a  function 
of  time  (Qrv(t)).  This  sketch  indicates  that  the  diastolic 
filling  time  period  (Td)  may  be  divided  into  three  phases. 
In  phase  1,  the  tricuspid  valve  opens  initiating  rapid  pas¬ 
sive  filling  from  the  right  atrium  to  the  right  ventricle. 
The  dynamics  involved  in  this  phase  are  principally  de¬ 
termined  by  the  right  atrial  and  ventricular  diastolic  com¬ 
pliances  and  the  small  tricuspid  resistance.  Phase  2  then 
follows  with  slower  passive  filling  from  the  systemic  veins 
to  the  right  ventricle  via  the  right  atrium  with  the  dy¬ 


namics  largely  determined  by  the  right  atrial  and  ven¬ 
tricular  diastolic  compliances,  the  tricuspid  resistance, 
and  the  relatively  larger  systemic  venous  resistance.  Fi¬ 
nally,  in  phase  3,  rapid  active  filling  from  the  contract¬ 
ing  right  atrium  to  the  right  ventricle  occurs  with  the 
dynamics  predominantly  determined  by  right  ventricular 
diastolic  compliance,  the  time-varying  compliance  of  the 
right  atrium  during  systole,  and  the  tricuspid  resistance. 

The  sketch  in  Fig.  1  indicates  that  Qrv(t)  rises  from  its 
end-systolic  value  (Q®®)  to  its  end-diastolic  value  (Q®(J) 
over  Td-  If  Td  were  infinitely  long,  Qrv(t)  would  rise  to  a 
maximum  value  (Q™ax)  that  is  usually  established  by  the 
product  of  the  right  ventricular  diastolic  compliance  and 
the  pressure  in  the  systemic  veins  which  is  essentially  con¬ 
stant  due  to  its  large  compliance.  A  unique  exponential 
function  may  be  fit  to  Td,  ,  Q®(J,  and  Q™ax  as  shown 
in  the  figure.  The  time  constant  of  this  exponential  is 
a  simple,  scalar  parameter  which  summarizes  ventricu¬ 
lar  diastolic  filling  dynamics.  We  refer  to  this  parameter 
as  the  effective  ventricular  diastolic  filling  time  constant 
(td).  This  time  constant  may  be  thought  of  as  the  prod¬ 
uct  of  a  lumped  resistance  and  a  lumped  compliance.  For 
example,  td  would  decrease  with  right  ventricular  dias¬ 
tolic  compliance  and  increase  with  tricuspid  resistance. 
Therefore,  td  could  be  utilized  as  a  means  for  monitor¬ 
ing  changes  in  these  resistances  and  compliances. 


Fig.  1.  Definition  of  a  scalar  parameter  td  which  summarizes 
ventricular  diastolic  filling  dynamics. 
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In  practice,  it  is  not  feasible  to  measure  directly  td 
because  of  the  difficulty  in  establishing  Q“ax.  However, 
in  the  next  section,  we  describe  one  method  to  measure 
td  indirectly  from  measurements  that  may  be  obtained 
noninvasively  in  humans. 


II.  Method 

The  method  is  based  on  the  assumption  that  a  small 
step  increase  in  HR  would  decrease  SV  in  the  steady-state 
due  solely  to  the  reduction  in  Td-  The  other  determi¬ 
nants  of  steady-state  SV,  namely  Q™ax  and  Q®®,  would 
be  relatively  unaffected  primarily  because  systemic  ve¬ 
nous  compliance  is  large  and  pulmonary  arterial  resis¬ 
tance  and  right  ventricular  end-systolic  compliance  are 
small.  By  further  assuming  a  single  time  constant  model 
of  ventricular  diastolic  filling  dynamics,  steady-state  SV 
may  then  be  written  according  to  the  following  equation: 

SV=SVmax{  (1) 


where  SVmax  =  Q™ax  —  Q®®  and  is  assumed  to  be  con¬ 
stant  for  small  HR  changes.  This  equation  indicates  that 
td  may  be  determined  from  two  steady-state  pairs  of  Td 
and  SV  measurements  through  solution  of  the  following 
nonlinear  equation: 


SVi  1  -  e-TD,i/TD 

SV2  _  1  -  e-TD,2/TD  ■ 


(2) 


The  two  steady-state  pairs  of  measurements  may  be 
obtained  by  slightly  varying  the  mean  HR  with  a  pacing 
protocol.  However,  they  may  also  be  determined  from 
a  single  pair  of  beat-to-beat  measurements  of  HR  and 
left  ventricular  SV  obtained  during  normal,  resting  con¬ 
ditions  in  which  the  fluctuations  in  these  variables  are 
small.  The  idea  here  is  to  establish  the  first  pair  with  the 
mean  values  of  HR  and  SV  (SVi  =  SV,  HRi  =  HR)  and 
the  second  pair  by  identifying  the  transfer  function  be¬ 
tween  the  fluctuations  in  HR  and  the  fluctuations  in  SV 
(HR— >SV)  and  computing  its  DC  value  (G)  -  the  steady- 
state  SV  response  to  a  unit  step  increase  in  HR.  That 
is,  SV2  =  SV  +  cthrG  and  HR2  =  HR  +  cthr,  where  cthr 
is  the  standard  deviation  of  HR.  For  small  changes  in 
HR,  Td,(i,2)  may  be  approximated  from  HR(1>2)  with  the 

scale  factor  (1  —  0.4/y'^^)  [2].  The  parameter  td  may 
then  be  determined  with  the  bisection  method  from  these 
two  steady-state  pairs  and  (2).  The  standard  deviation 
of  td  may  be  analogously  estimated  with  the  standard 
deviation  of  the  estimate  of  G. 

In  order  to  identify  reliably  HR— >SV,  it  is  necessary 
to  measure  and  account  for  ABP  fluctuations  which  also 
contribute  to  left,  ventricular  SV  fluctuations  through 
mechanical  afterload  and  baroreflex  effects.  This  may 
be  achieved  by  simultaneously  identifying  the  transfer 
function  between  ABP  fluctuations  and  SV  fluctuations 
(ABP— >SV)  from  the  beat-to-beat  measurements  as  de¬ 
picted  in  the  block  diagram  in  Fig.  2.  The  term  Nsv 


HR— ►SV 

4- 


ABP— SV 

Fig.  2.  Identification  method  for  establishing  the  two  steady-state 
pairs  of  HR  and  SV  values. 


in  the  block  diagram  is  also  estimated  from  the  beat-to- 
beat  measurements  and  represents  the  residual  variability 
in  SV  not  accounted  for  by  HR  or  ABP  fluctuations. 

The  block  diagram  in  Fig.  2  is  mathematically  rep¬ 
resented  by  an  autoregressive  moving  average  (ARMA) 
difference  equation  of  the  following  form: 

m  n 

sv(t)  =  5>  SV{t  -i)  +  ^2  biHR(t  -  i) 

2=1  2=0 

P 

+  CiABP(t  -  i)  +  Wsv(t),  (3) 

i—0 

where  t  is  discrete  time  and  Wsv  is  the  residual  er¬ 
ror.  The  three  sets  of  parameters  {a,,  bp  c,}  completely 
specify  the  transfer  functions  characterizing  HR— >SV  and 
ABP— »SV,  while  the  residual  error  together  with  the  set 
of  parameters  {a;}  fully  define  Nsv  [3].  The  values  of  the 
parameters  are  determined  from  six  minute  segments  of 
zero-mean  HR,  ABP,  and  SV  sampled  at  0.5  Hz  by  min¬ 
imizing  the  variance  of  the  residual  error  in  conjunction 
with  an  ARMA  parameter  reduction  algorithm  [3].  The 
mean  of  the  estimate  of  G  is  then  given  by  1  Vjvff  b‘a  , 
while  its  standard  deviation  may  be  approximated  as  de¬ 
scribed  in  [2],  [3]. 

III.  Cardiovascular  system  model 

A  computational  model  of  the  cardiovascular  system 
provides  a  convenient  means  for  evaluating  the  method, 
because  it  would  be  possible  to  establish  the  actual  td 
according  to  Fig.  1.  To  this  end,  we  constructed  a  compu¬ 
tational  model  of  the  human  cardiovascular  system  that 
is  capable  of  reasonably  emulating  experimental  beat-to- 
beat  hemodynamic  variability.  The  model  includes  three 
major  components:  a  heart  and  circulation,  a  short-term 
regulatory  system,  and  resting  physiologic  perturbations. 

The  heart  and  circulation  is  a  lumped  parameter  model 
which  includes  eight  compartments  representing  the  left 
atrium  and  ventricle,  the  right  atrium  and  ventricle,  the 
systemic  arteries  and  veins,  and  the  pulmonary  arter- 
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ies  and  veins.  Each  compartment  consists  of  a  linear  or 
nonlinear  resistance  and  a  linear  or  nonlinear  compliance 
with  an  associated  dead  volume.  The  compliances  of  the 
ventricular  and  atrial  compartments  vary  over  each  car¬ 
diac  cycle  so  as  to  drive  blood  flow  in  the  model. 

The  short-term  regulatory  system  consists  of  arterial 
and  cardiopulmonary  baroreflexes  as  well  as  a  direct  neu¬ 
ral  coupling  mechanism  between  respiration  and  heart 
rate.  The  arterial  baroreflex  maintains  ABP  through  the 
control  of  HR,  ventricular  end-systolic  compliances,  sys¬ 
temic  arterial  resistance,  and  systemic  venous  dead  vol¬ 
ume,  while  the  cardiopulmonary  baroreflex  senses  right 
atrial  transmural  pressure  and  manipulates  systemic  ar¬ 
terial  resistance  and  systemic  venous  dead  volume. 

The  resting  physiologic  perturbations  include  fixed- 
rate  respiratory  activity  which  impinges  upon  the  heart 
and  circulation  through  intrathoracic  pressure  -  the  ref¬ 
erence  pressure  of  the  pulmonary,  atrial,  and  ventricu¬ 
lar  compartments  -  as  well  as  the  direct  neural  coupling 
mechanism.  The  resting  physiologic  perturbations  also 
include  a  broadband  stochastic  disturbance  to  systemic 
arterial  resistance  which  is  primarily  responsible  for  gen¬ 
erating  low  frequency  (0.04-0.1  Hz)  fluctuations  through 
a  model  resonance  phenomenon  as  well  as  a  1/f  stochas¬ 
tic  disturbance  to  HR  which  initiates  very  low  frequency 
(<0.04  Hz)  variability.  Further  details  of  the  computa¬ 
tional  model  are  provided  in  [2]. 

IV.  Model  results 

Fig.  3  illustrates  the  results  of  the  application  of  the 
method  to  the  computational  model  generated  data  in 
terms  of  the  estimated  versus  actual  td-  The  actual  td 
was  varied  by  adjusting  the  right  ventricular  diastolic 
compliance  and  the  tricuspid  resistance  of  the  computa¬ 
tional  model.  The  estimates  are  shown  in  terms  of  their 
meanistandard  deviation  as  determined  from  20  differ¬ 
ent  realizations  of  computational  model  generated  data. 
The  solid  trace  is  the  identity  line,  and  the  dashed  trace 
is  the  line  of  best  fit.  The  results  here  suggest  that  td 
may  be  reliably  estimated  from  experimental  data  to  the 
extent  that  the  computational  model  coincides  with  real¬ 
ity. 

Fig.  4  illustrates  an  example  of  how  the  method  would 
perform  in  the  case  of  cardiac  tamponade  which  was 
implemented  in  the  computational  model  by  decreasing 
both  ventricular  diastolic  compliances.  The  figure  here 
illustrates  that  the  estimated  rD  (mean±  standard  de¬ 
viation)  successfully  decreases  with  right  ventricular  di¬ 
astolic  compliance  (C°).  Fig.  5  illustrates  an  example 
of  how  the  method  would  perform  in  the  case  of  mi¬ 
tral  stenosis  which  was  implemented  in  the  computational 
model  by  increasing  the  mitral  resistance  (Rmitrai)-  The 
figure  illustrates  that  the  estimated  td  (meanistandard 
deviation)  successfully  increases  with  Rmitrai-  Impor¬ 
tantly,  in  order  to  monitor  changes  in  and  Rmit  by 
conventional  methods,  invasive  measurements  would  be 


Fig.  3.  Computational  model  results.  See  text  for  details. 


required. 


V.  Method  limitations  and  applications 

Although  the  method  performs  quite  well  with  respect 
to  the  computational  model  generated  data,  its  potential 
applications  are  somewhat  limited.  Firstly,  the  definition 
of  td  as  depicted  in  Fig.  1  is  not  applicable  for  sufficiently 
fast  ventricular  diastolic  filling  dynamics  such  that  the 
atrial  contraction  causes  Q®((  to  exceed  Q“ax.  In  this 
case,  the  estimated  change  in  td  from  one  measurement 
period  to  the  next  may  be  misleading.  However,  the  ac¬ 
companying  change  in  the  SV  measurement  may  aid  in 
clarifying  the  actual  change  in  td- 
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Fig.  4.  Computational  model  example  of  cardiac  tamponade  mon¬ 
itoring.  See  text  for  details. 
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VI.  Conclusions 

In  this  paper,  we  presented  a  method  to  estimate  a 
parameter  summarizing  the  diastolic  filling  dynamics  of 
the  ventricle  which  governs  the  effects  of  Td  on  SV.  The 
method  requires  beat-to-beat  measurements  of  HR,  SV, 
and  ABP,  each  of  which  may  be  obtained  noninvasively 
in  humans.  We  validated  the  method  against  data  gener¬ 
ated  from  a  computational  model  of  the  human  cardiovas¬ 
cular  system.  The  computational  model-based  validation 
motivates  experimental  evaluation  of  the  method.  This 
may  be  achieved,  for  example,  by  reducing  the  ventricu¬ 
lar  diastolic  compliances  in  an  animal  model  of  cardiac 
tamponade  and  assessing  whether  the  td  estimate  is  able 
to  track  these  reductions. 
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Fig.  5.  Computational  model  example  of  mitral  stenosis  monitor¬ 
ing.  See  text  for  details. 


Secondly,  the  method  could  not  be  used  for  diagnostic 
purposes.  That  is,  a  change  in  td  could  be  due  to  a 
change  in  the  state  of  either  a  resistance  or  compliance. 
A  tq  change  could  also  be  due  to  a  shift  in  operating  point 
of  either  a  resistance  or  compliance,  since  td  characterizes 
ventricular  diastolic  filling  dynamics  over  a  narrow  range 
of  operating  points  established  by  the  size  of  the  analyzed 
hemodynamic  fluctuations.  However,  once  the  diagnosis 
is  made,  td  could  then  be  used  for  monitoring  purposes. 

Finally,  the  estimated  td  actually  represents  the  dias¬ 
tolic  filling  dynamics  of  the  ventricle  which  govern  the 
effects  of  Td  on  SV.  Assuming  that  the  right  ventricu¬ 
lar  diastolic  compliance  and  the  lumped  inflow  resistance 
are  both  larger  than  their  left  ventricular  counterparts 
(as  is  the  case  in  the  computational  model),  then  td  nor¬ 
mally  characterizes  right  ventricular  diastolic  filling  dy¬ 
namics.  This  implies  that  the  method  could  not  be  used 
to  monitor  isolated  left  ventricular  diastolic  compliance 
reductions  ( e.g .,  hypertrophic  cardiomyopathy,  hypertro¬ 
phy  due  to  hypertension,  left  ventricular  ischemia),  be¬ 
cause  the  slower  right  ventricle  would  govern  the  effects 
of  Td  on  SV  in  this  case.  However,  the  method  could 
be  used  to  monitor  simultaneous  reductions  in  both  ven¬ 
tricular  diastolic  compliances  (e.g.,  cardiac  tamponade, 
restrictive  cardiomyopathy),  isolated  right  ventricular  di¬ 
astolic  compliance  reductions  (e.g.  hypertrophy  due  pul¬ 
monary  hypertension),  and  lumped  inflow  resistance  in¬ 
creases  (e.g.,  stenosis  of  the  tricsupid  and  mitral  valves). 
We  acknowledge  that  it  may  be  possible  to  character¬ 
ize  td  and  its  associated  operating  point  for  each  ventri¬ 
cle  with  beat-to-beat  ventricular  volume  measurements 
obtained  with  B-mode  echocardiography.  However,  this 
approach  also  has  its  disadvantages  including  exorbitant 
image  processing. 
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